% wrote by deepseek-R1 & TYY
function dX = ode_system(X_prev, U_prev, h_t_prev, grd_t_prev, params)
    % 解包参数
    abc = params.abc; m = params.m; g = params.g; y = params.y;
    cs1 = params.cs1; cs2 = params.cs2; cs3 = params.cs3; cs4 = params.cs4; cs5 = params.cs5;
    cr1 = params.cr1; cr2 = params.cr2; cr3 = params.cr3; cr4 = params.cr4; cr5 = params.cr5;
    Cs = params.Cs; Cr = params.Cr; d1 = params.d1; d2 = params.d2; k5 = params.k5; k7 = params.k7;
    mt_num = params.mt_num;

    % 速度导数 (动力学方程)
    ddv = ([1 -1]*U_prev - (abc(1)+abc(2)*X_prev(1)*3.6 + abc(3)*(X_prev(1)^2)*3.6^2)*m*g/1000 ...
           - grd_t_prev*m*g/1000) / (m*(1+y));
    
    % 位移导数 (速度)
    dds = X_prev(1);  % 修正原代码中的错误：位移导数为速度
    
    % 定子温度导数
    ddtheta_s = cs1*X_prev(3) + cs2*X_prev(4) + cs3*(X_prev(1)^2) + ...
                cs4*(U_prev(1)^2)/mt_num^2 + cs5*(20 -6*h_t_prev/1000) - ...
                (1/Cs)*(d1 + k5*h_t_prev);
    
    % 转子温度导数
    ddtheta_r = cr1*X_prev(3) + cr2*X_prev(4) + cr3*(X_prev(1)^2) + ...
                cr4*(U_prev(1)^2)/mt_num^2 + cr5*(20 -6*h_t_prev/1000) - ...
                (1/Cr)*(d2 + k7*h_t_prev);
    
    dX = [ddv; dds; ddtheta_s; ddtheta_r];
end